# ------------------------------------- #
# Replication code for:
#
# Pomeroy, Caleb "Hawks Become Us: The Sense of Power and Militant Foreign Policy Attitudes," Security Studies.
#
# This script reproduces the key results from the paper.
# ------------------------------------- #

# --- set your working directory --- #
setwd("~/Downloads/hawks_replication/")

# --- libraries --- #
library(readr)
library(ggplot2)
library(texreg)
library(mediation)
library(plyr)
library(scales)
library(caret)
library(ggeffects)
library(ggh4x)


# --- import data --- #
# --- app data
app_survey <- read_csv("data/app_survey.csv")
app_survey$gender <- factor(app_survey$gender, levels = c("non_male", "male"))
# NB: to ease interpretation, sense of power (higher = greater feeling of power) and ideology (higher = more conservative) are scaled to mean 0 and SD 1; 
# MI is scaled to 0-1; age is retained on a numeric scale. 

# --- russia data
russia_survey <- read_csv("data/russia_survey.csv")
russia_survey$gender <- factor(russia_survey$gender, levels = c("female", "male"))
russia_survey$education <- factor(russia_survey$education, levels = c("bachelors", "less_bachelors"))
# NB: each of the following scaled to mean 0 and SD 1: sense of power (higher = greater feeling of power), national attachment, harm/care, fairness/reciprocity, ingroup/loyalty, authority/respect;
# MI is scaled to 0-1; age is retained on a numeric scale; the maritime security DV ("boat_threat") is retained on a Likert scale (higher = more unilaterally aggressive)

# --- qualtrics US data (survey #1, which included in-depth scenario)
qualtrics_survey_scenario <- read_csv("data/us_survey1.csv")
qualtrics_survey_scenario$gender <- factor(qualtrics_survey_scenario$gender, levels = c("non_male", "male"))
# NB: each of the following scaled to mean 0 and SD 1: sense of power (higher = greater feeling of power), national attachment, political ideology (higher = more conservative); for party ID ("pid"), "non_republican" captures moderates and democrats;
# MI, preventive strikes ("iran_strike"), and nuke support ("iran_nukes") are scaled to 0-1; age was collected and retained in bins (1 = 18-24, 2 = 25-34, 3 = 35-44, 4 = 45-54, 5 = 55-64, 6 = 65-74, 7 = 75-84, 8 = 85+)

# --- china data
china_survey <- read_csv("data/china_survey.csv")
china_survey$gender <- factor(china_survey$gender, levels = c("female", "male"))
# NB: each of the following scaled to mean 0 and SD 1: sense of power (higher = greater feeling of power), national attachment, harm/care, fairness/reciprocity, ingroup/loyalty, authority/respect;
# binding foundations is simply the sum of authority/respect and ingroup/loyalty scores
# MI, aggressive intentions ("aggressive_intentions"), and US threat perception ("us_threats") are scaled to 0-1; age was collected and retained in bins (1 = 18-24, 2 = 25-34, 3 = 35-44, 4 = 45-54, 5 = 55-64, 6 = 65-74, 7 = 75-84, 8 = 85+)

# --- qualtrics US data (survey #2, which included multiple measures of individual differences)
qualtrics_survey_inddiffs <- read_csv("data/us_survey2.csv")
qualtrics_survey_inddiffs$gender <- factor(qualtrics_survey_inddiffs$gender, c("non_male", "male"))
qualtrics_survey_inddiffs$education <- factor(qualtrics_survey_inddiffs$education, levels = c("bachelors", "less_bachelors"))
# NB: each of the following scaled to mean 0 and SD 1: sense of power (higher = greater feeling of power), national attachment, harm/care, fairness/reciprocity, ingroup/loyalty, authority/respect,
# extraversion, agreeableness, conscientiousness, neuroticism, openness, rwa, sdo;
# MI is scaled to 0-1; age was collected and retained in bins (1 = 18-24, 2 = 25-34, 3 = 35-44, 4 = 45-54, 5 = 55-64, 6 = 65-74, 7 = 75-84, 8 = 85+)
# the individual "sense_power_" and "mi_" items are retained on original Likert scales, as described in the paper

# --- prolific US data (survey experiment)
prolific_survey <- read_csv("data/prolific_survey.csv")
prolific_survey$gender <- factor(prolific_survey$gender, levels = c("non_male", "male"))
prolific_survey$education <- factor(prolific_survey$education, levels = c("less_bachelors", "bachelors"))
# NB: each of the following represented with unidimensional factor scores: MI ("mi_factor"), sense of power ("sense_power_factor")
# each of the following scaled to mean 0 and SD 1: ideology, pid, national attachment.
# age was collected and retained in bins (1 = 18-24, 2 = 25-34, 3 = 35-44, 4 = 45-54, 5 = 55-64, 6 = 65-74, 7 = 75-84)



# --- Sense of State Power Explains Militant Internationalism --- #
# These are the main models from the paper, showing that sense of power explains positive variance in MI

summary(mi_us_simple <- lm(mil_int ~ sense_power, data = qualtrics_survey_inddiffs))
summary(mi_us <- lm(mil_int ~ sense_power + nat_attach + education + gender + age + 
                      harm_care + fairness_recip + ingroup_loyalty + authority_respect +
                      extraversion + agreeableness + conscientiousness + neuroticism + openness +
                      rwa + sdo, data = qualtrics_survey_inddiffs))

summary(mi_china_simple <- lm(mi_war ~ sense_power, data = china_survey))
summary(mi_china <- lm(mi_war ~ sense_power + nat_attach +  education + gender + age +
                         harm_care + fairness_recip + ingroup_loyalty + authority_respect, data = china_survey))

summary(mi_russia_simple <- lm(mil_int ~ sense_power, data = russia_survey))
summary(mi_russia <- lm(mil_int ~ sense_power + nat_attach + education + gender + age +
                          harm_care + fairness_recip + ingroup_loyalty + authority_respect, data = russia_survey))

# --- Table A4:
screenreg(list(mi_us_simple, mi_us, mi_china_simple, mi_china, mi_russia_simple, mi_russia),
          booktabs = T,
          caption.above = T,
          caption = "OLS Results: Sense of Power Explains Militant Orientation",
          custom.coef.names = c("(Intercept)",
                                "Sense of State Power",
                                "National Attachment",
                                "Less than College",
                                "Male",
                                "Age",
                                "Harm/Care",
                                "Fairness/Reciprocity",
                                "Ingroup/Loyalty",
                                "Authority/Respect",
                                "Extraversion",
                                "Agreeableness",
                                "Conscientiousness",
                                "Neuroticism",
                                "Openness",
                                "Right-Wing Authoritarianism",
                                "Social Dominance Orientation"),
          stars = c(.001, .01, .05, .10),
          symbol = "\\wedge",
          custom.model.names = c("(1)",
                                 "(2)",
                                 "(3)",
                                 "(4)",
                                 "(5)",
                                 "(6)"))

mod_list <- list("us_public" = mi_us,
                 "chinese_public" = mi_china,
                 "russian_public" = mi_russia)
results_df <- data.frame()
for(i in 1:length(mod_list)){
  mod_i <- mod_list[[i]]
  results_df <- rbind(results_df,
                      data.frame(ci_95_high = confint(mod_i, level = .95)[,2],
                                 ci_95_low = confint(mod_i, level = .95)[,1],
                                 ci_90_high = confint(mod_i, level = .90)[,2],
                                 ci_90_low = confint(mod_i, level = .90)[,1],
                                 est = coef(mod_i),
                                 variable = names(coef(mod_i)),
                                 model = names(mod_list)[i]))
  
}

results_df$variable <- revalue(results_df$variable, c("sense_power" = "Sense of Power",
                                                      "nat_attach" = "National Attachment",
                                                      "ingroup_loyalty" = "Ingroup/Loyalty",
                                                      "harm_care" = "Harm/Care",
                                                      "fairness_recip" = "Fairness/Reciprocity",
                                                      "authority_respect" = "Authority/Respect",
                                                      "gendermale" = "Male",
                                                      "rwa" = "Right-Wing Authoritarianism",
                                                      "sdo" = "Social Dominance Orientation",
                                                      "age" = "Age",
                                                      "extraversion" = "Extraversion",
                                                      "agreeableness" = "Agreeableness",
                                                      "conscientiousness" = "Conscientiousness", 
                                                      "neuroticism" = "Neuroticism",
                                                      "openness" = "Openness",
                                                      "educationless_bachelors" = "No College")) 

results_df$model <- factor(results_df$model, levels = c("us_public", "chinese_public", "russian_public"))
levels(results_df$model) <- c("American Public", "Chinese Public", "Russian Public")
results_df <- results_df[!results_df$variable%in%c("(Intercept)", "Age", "Male", "No College", "National Attachment"),]
results_df$variable <- 
  factor(results_df$variable, rev(c("Sense of Power",
                                    "Ingroup/Loyalty", "Authority/Respect", "Harm/Care", "Fairness/Reciprocity", 
                                    "Extraversion", "Agreeableness", "Conscientiousness", "Neuroticism", "Openness",
                                    "Right-Wing Authoritarianism", "Social Dominance Orientation")))
plot_col <- ifelse(results_df$variable == "Sense of Power", "#54545d", "#99998e")

# --- Figure 1:
ggplot(results_df, aes(x = est, y = variable, color = plot_col)) +
  geom_vline(xintercept = 0, linetype = "solid", color = "black", size = .8) +
  geom_linerange(aes(xmin = ci_95_low, xmax = ci_95_high), size = 1, color = "black") +
  geom_linerange(aes(xmin = ci_90_low, xmax = ci_90_high), size = 2.8) +
  geom_point(size = 5, color = "black") +
  geom_point(size = 4) +
  facet_wrap(~model, nrow = 1) +
  theme_bw() +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(x = "Estimate", y = NULL) +
  scale_color_manual(values = plot_col) + 
  theme(axis.title.x = element_text(size = 18, margin = margin(t=15)),
        axis.text.y = element_text(size = 16),
        axis.text.x = element_text(size = 13),
        strip.text = element_text(size = 14, color = "gray94"),
        legend.position = "none",
        panel.spacing = unit(1.2, "lines"),
        panel.background = element_rect(fill = "#f6f6f5"),
        strip.background = element_rect(fill = "#54544c", color = "black", size = 1.5),
        plot.margin = margin(r=15, l=15, t= 15))


# --- Sense of State Power Explains Discrete Attitudes in Country-Specific Scenarios --- #
# Beyond the primary MI DV above, three surveys included in-depth country-specific security scenarios to assess whether 
# the sense of power explains positive variance in discrete attitudes beyond generalized hawkishness.

# --- The US public and Iran's nuclear program 
summary(iran_power_strike <- lm(iran_strike ~ sense_power + ideology + gender + age + white + pid + nat_attach, data = qualtrics_survey_scenario))
summary(iran_power_nukes <- lm(iran_nukes ~ sense_power + ideology + gender + age + white + pid + nat_attach, data = qualtrics_survey_scenario))

mod_list <- list("iran_power_strike" = iran_power_strike,
                 "iran_power_nukes" = iran_power_nukes)

# --- Table A12
screenreg(list(iran_power_strike, iran_power_nukes),
       booktabs = T,
       caption.above = T,
       caption = "OLS Results: Sense of American Power Explains Discrete Security Attitudes",
       custom.coef.names = c("(Intercept)",
                             "Sense of American Power",
                             "Conservative",
                             "Male",
                             "Age",
                             "White",
                             "Republican",
                             "National Attachment"),
       stars = c(.001, .01, .05, .10),
       symbol = "\\wedge",
       label = "table:ols_iran_scenario",
       custom.model.names = c("Risky Preventive Strike",
                              "Nuclear Weapons Use"))

results_df_us <- data.frame()
for(i in 1:length(mod_list)){
  mod_i <- mod_list[[i]]
  results_df_us <- rbind(results_df_us,
                      data.frame(ci_95_high = confint(mod_i, level = .95)[,2],
                                 ci_95_low = confint(mod_i, level = .95)[,1],
                                 ci_90_high = confint(mod_i, level = .90)[,2],
                                 ci_90_low = confint(mod_i, level = .90)[,1],
                                 est = coef(mod_i),
                                 variable = names(coef(mod_i)),
                                 model = names(mod_list)[i]))
}
results_df_us$variable <- revalue(results_df_us$variable, c("sense_power" = "Sense of Power",
                                                      "ideology" = "Conservative",
                                                      "gendermale" = "Male",
                                                      "age" = "Age",
                                                      "whitewhite" = "White",
                                                      "pidrepublican" = "Republican",
                                                      "nat_attach" = "National Attachment")) 

# we'll wrap the three surveys' figures together below, just setting up formatting for now:
results_df_us$model <- factor(results_df_us$model, levels = c("iran_power_strike", "iran_power_nukes"))
levels(results_df_us$model) <- c("Support for Risky\nPreventive Strike", "Support for Nuclear\nWeapons Use")
results_df_us <- results_df_us[!results_df_us$variable%in%c("(Intercept)", "National Attachment", "White", "Republican", "Conservative"),]
results_df_us$variable <- factor(results_df_us$variable, rev(c("Sense of Power", "Cooperative Internationalism", "Male", "Age")))
results_df_us$plot_col <- ifelse(results_df_us$variable == "Sense of Power", "#54545d", "#99998e")
results_df_us$country <- "us"


# --- The Chinese public and competition in the South China Sea 
summary(china_power_intentions <- lm(aggressive_intentions ~ sense_power + nat_attach +  education + gender + age +
                                       harm_care + fairness_recip + ingroup_loyalty + authority_respect, data = china_survey))
summary(china_power_threats <- lm(us_threats ~ sense_power + nat_attach +  education + gender + age +
                                    harm_care + fairness_recip + ingroup_loyalty + authority_respect, data = china_survey))

# --- Table A13:
screenreg(list(china_power_intentions, china_power_threats),
       booktabs = T,
       caption.above = T,
       caption = "OLS Results: Sense of Chinese Power Explains Discrete Security Attitudes",
       custom.coef.names = c("(Intercept)",
                             "Sense of Chinese Power",
                             "National Attachment",
                             "No College Degree",
                             "Male",
                             "Age",
                             "Harm/Care",
                             "Fairness/Reciprocity",
                             "Ingroup/Loyalty",
                             "Authority/Respect"),
       stars = c(.001, .01, .05, .10),
       symbol = "\\wedge",
       label = "table:ols_china_scenario",
       custom.model.names = c("Countries Harbor Aggressive Intentions",
                              "U.S. Threatens China's Security"))

mod_list <- list("china_power_threats" = china_power_threats,
                 "china_power_intentions" = china_power_intentions)

results_df_china <- data.frame()
for(i in 1:length(mod_list)){
  mod_i <- mod_list[[i]]
  results_df_china <- rbind(results_df_china,
                      data.frame(ci_95_high = confint(mod_i, level = .95)[,2],
                                 ci_95_low = confint(mod_i, level = .95)[,1],
                                 ci_90_high = confint(mod_i, level = .90)[,2],
                                 ci_90_low = confint(mod_i, level = .90)[,1],
                                 est = coef(mod_i),
                                 variable = names(coef(mod_i)),
                                 model = names(mod_list)[i]))
}
results_df_china$variable <- revalue(results_df_china$variable, c("sense_power" = "Sense of Power",
                                                      "educationless_bachelors" = "No College Degree",
                                                      "gendermale" = "Male",
                                                      "age" = "Age",
                                                      "harm_care" = "Harm/Care",
                                                      "fairness_recip" = "Fairness/Reciprocity",
                                                      "ingroup_loyalty" = "Ingroup/Loyalty",
                                                      "authority_respect" = "Authority/Respect",
                                                      "nat_attach" = "National Attachment")) 

results_df_china$model <- factor(results_df_china$model, levels = c("china_power_intentions", "china_power_threats"))
levels(results_df_china$model) <- c("Countries Harbor\nAggressive Intentions", "U.S. Threatens\nChina's Security")
results_df_china <- results_df_china[!results_df_china$variable%in%c("(Intercept)", "Important to Perspective-Take", "National Attachment", "No College Degree", "Harm/Care", "Fairness/Reciprocity"),]
results_df_china$variable <- factor(results_df_china$variable, rev(c("Sense of Power", "Ingroup/Loyalty", "Authority/Respect", "Male", "Age")))
results_df_china$plot_col <- ifelse(results_df_china$variable == "Sense of Power", "#54545d", "#99998e")
results_df_china$country <- "china"

# --- The Russian public and maritime security
russia_survey$boat_threat_binary <- ifelse(russia_survey$boat_threat == 5, 1, 0)
summary(boat_russia <- glm(boat_threat_binary ~ sense_power + nat_attach + education + gender + age +
                             harm_care + fairness_recip + ingroup_loyalty + authority_respect, 
                           data = russia_survey, family = "binomial"))

# result is robust to retaining the DV on a numeric scale, too:
# summary(lm(boat_threat ~ sense_power + nat_attach + education + gender + age + harm_care + fairness_recip + ingroup_loyalty + authority_respect, data = russia_survey))

# --- Table A11:
screenreg(boat_russia,
       booktabs = T,
       caption.above = T,
       caption = "OLS Results: Sense of Russian Power Explains Discrete Security Attitudes",
       custom.coef.names = c("(Intercept)",
                             "Sense of Russian Power",
                             "National Attachment",
                             "No College Degree",
                             "Male",
                             "Age",
                             "Harm/Care",
                             "Fairness/Reciprocity",
                             "Ingroup/Loyalty",
                             "Authority/Respect"),
       stars = c(.001, .01, .05, .10),
       symbol = "\\wedge",
       label = "table:logit_russia_scenario",
       custom.model.names = c("Strike the Foreign Warship"))

mod_list <- list("boat_russia" = boat_russia)
results_df_russia <- data.frame()
for(i in 1:length(mod_list)){
  mod_i <- mod_list[[i]]
  results_df_russia <- rbind(results_df_russia,
                      data.frame(ci_95_high = confint(mod_i, level = .95)[,2],
                                 ci_95_low = confint(mod_i, level = .95)[,1],
                                 ci_90_high = confint(mod_i, level = .90)[,2],
                                 ci_90_low = confint(mod_i, level = .90)[,1],
                                 est = coef(mod_i),
                                 variable = names(coef(mod_i)),
                                 model = names(mod_list)[i]))
}
results_df_russia$variable <- revalue(results_df_russia$variable, c("sense_power" = "Sense of Power",
                                                      "educationless_bachelors" = "No College Degree",
                                                      "gendermale" = "Male",
                                                      "age" = "Age",
                                                      "harm_care" = "Harm/Care",
                                                      "fairness_recip" = "Fairness/Reciprocity",
                                                      "ingroup_loyalty" = "Ingroup/Loyalty",
                                                      "authority_respect" = "Authority/Respect",
                                                      "nat_attach" = "National Attachment")) 

results_df_russia$model <- factor(results_df_russia$model, levels = c("boat_russia"))
levels(results_df_russia$model) <- c("Strike the Foreign\nWarship")
results_df_russia <- results_df_russia[!results_df_russia$variable%in%c("(Intercept)", "Important to Perspective-Take", "National Attachment", "No College Degree", "Harm/Care", "Fairness/Reciprocity"),]
results_df_russia$variable <- factor(results_df_russia$variable, rev(c("Sense of Power", "Ingroup/Loyalty", "Authority/Respect", "Male", "Age")))
results_df_russia$plot_col <- ifelse(results_df_russia$variable == "Sense of Power", "#54545d", "#99998e")
results_df_russia$country <- "russia"

# combine the various figures
results_df <- rbind(results_df_us, results_df_china, results_df_russia)
results_df$country <- factor(results_df$country, levels = c("us", "china", "russia"))
levels(results_df$country) <- c("American Public", "Chinese Public", "Russian Public")
results_df$variable <- factor(results_df$variable, levels = c("Authority/Respect", "Ingroup/Loyalty", 
                                                              "Age", "Male", "Sense of Power"))
p <- 
  ggplot(results_df, aes(x = est, y = variable, color = plot_col)) +
  geom_vline(xintercept = 0, linetype = "solid", color = "black", size = .8) +
  geom_linerange(aes(xmin = ci_95_low, xmax = ci_95_high), size = 1, color = "black") +
  geom_linerange(aes(xmin = ci_90_low, xmax = ci_90_high), size = 2.8) +
  geom_point(size = 5, color = "black") +
  geom_point(size = 4) +
  facet_nested(~country+model, scales = "free_x") +
  theme_bw() +
  labs(x = "Estimate", y = NULL) +
  scale_x_continuous(n.breaks = 4) +
  scale_color_manual(values = plot_col) + 
  theme(axis.title.x = element_text(size = 17, margin = margin(t=15)),
        axis.text.y = element_text(size = 14),
        axis.text.x = element_text(size = 12),
        strip.text = element_text(size = 12, color = "gray94"),
        legend.position = "none",
        panel.spacing = unit(.9, "lines"),
        panel.background = element_rect(fill = "#f6f6f5"),
        strip.background = element_rect(fill = "#54544c", color = "black", size = 1.5),
        plot.margin = margin(r=75, l=55, t= 15))

# --- Figure 2:
# simply switches up the sub-panel background color:
{g <- ggplot_gtable(ggplot_build(p))
  strip_it <- which(grepl('strip-t', g$layout$name))
  colors <- c(rep("gray94", 3), rep("black", 5))
  fills <- c(rep("#54544c", 3),rep("gray95", 5))
  k <- 1
  for (i in strip_it) {
    j <- which(grepl("text", g$grobs[[i]]$grobs[[1]]$childrenOrder))
    g$grobs[[i]]$grobs[[1]]$children[[j]]$children[[1]]$gp$col <- colors[k]
    t <- which(grepl("rect", g$grobs[[i]]$grobs[[1]]$childrenOrder))
    g$grobs[[i]]$grobs[[1]]$children[[t]]$gp$fill <- fills[k]
    k <- k+1
  }
  grid::grid.draw(g)
  }



# --- Why Does Power Explain Hawkishness? Moderation Analyses--- #
# The above results suggest that the sense of power explains hawkishness, measured as both generalized MI and 
# discrete, scenario-specific attitudes. Who drives these results? That is, does power reveal, corrupt, or both?

# --- US moderation analysis: power reveals conservatives, corrupts liberals (MI)
app_survey$ideology_factor <- ifelse(app_survey$ideology <= median(app_survey$ideology), "liberal", "conservative")
app_survey$ideology_factor <- factor(app_survey$ideology_factor, levels = c("conservative", "liberal"))

summary(app_moderation <- lm(mi_war ~ sense_power*ideology_factor + gender + age + urban_rural + education + white, data = app_survey))
app_moderation <- ggpredict(app_moderation, c("sense_power", "ideology_factor")) 

plot_df_us_ideology_mi <- 
  rbind(data.frame(mil_int = app_moderation$predicted,
                   iv = app_moderation$group,
                   power = app_moderation$x,
                   ci_high = app_moderation$conf.high,
                   ci_low = app_moderation$conf.low))

levels(plot_df_us_ideology_mi$iv ) <- c("Conservative", "Liberal")
plot_df_us_ideology_mi$country <- "American Public"
# becomes panel A of Figure 3, to be continued below
plot_us_ideology_mi <- 
  ggplot(plot_df_us_ideology_mi, aes(x = power, y = mil_int, group = iv)) +
  geom_ribbon(aes(ymax = ci_high, ymin = ci_low), alpha = .45, fill = "gray70") +
  geom_line(aes(color = iv), size = 2) +
  labs(y = "Militant Internationalism", x = "Sense of American Power", colour = "Political Ideology") +
  scale_color_manual(values = c("#99b0b0", "#54545d")) +
  theme_bw() +
  facet_wrap(~country) +
  guides( colour = guide_legend(reverse = T)) +
  theme(axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14),
        axis.title.y = element_text(size = 16, margin = margin(r=10)),
        axis.title.x = element_text(size = 16, margin = margin(t=10)),
        legend.text = element_text(size = 13),
        legend.title = element_text(size = 13),
        strip.background = element_rect(fill = "#54544c", color = "black", size = 1.5),
        panel.background = element_rect(fill = "#f6f6f5"),
        panel.spacing = unit(1.3, "lines"),
        strip.text = element_text(size = 14, color = "gray94"),
        plot.margin = margin(20,20,20,20))


# --- Table A7:
# The above ideology results are dichotomized for ease of presentation in the main text;
# these appendix models represent ideology using the original scaled values
summary(app_moderation_simple <- lm(mi_war ~ sense_power*ideology, data = app_survey))
summary(app_moderation <- lm(mi_war ~ sense_power*ideology + gender + age + urban_rural + education + white, data = app_survey))

screenreg(list(app_moderation_simple, app_moderation),
       booktabs = T,
       caption.above = T,
       caption = "OLS Results: Sense of Power Reveals Conservatives, Corrupts Liberals",
       custom.coef.names = c("(Intercept)",
                             "Sense of State Power",
                             "Conservative",
                             "Sense of State Power $\\times$ Conservative",
                             "Male",
                             "Age",
                             "Urban Resident",
                             "Less than College",
                             "White"),
       stars = c(.001, .01, .05, .10),
       symbol = "\\wedge",
       custom.model.names = c("(1)",
                              "(2)"))

# --- for US
qualtrics_survey_scenario$ideology_liberal <- qualtrics_survey_scenario$ideology*-1 # multiply conservative scores by -1 to interpret as liberal scores (i.e. less conservative)
summary(us_power_strike_numeric <- lm(iran_strike ~ sense_power*ideology_liberal + ideology_liberal + gender + age + white + pid + nat_attach,
                                      qualtrics_survey_scenario))
summary(us_power_nukes_numeric <- lm(iran_nukes ~ sense_power*ideology_liberal + ideology_liberal + gender + age + white + pid + nat_attach,
                                     qualtrics_survey_scenario))

# --- Table A8
screenreg(list(us_power_strike_numeric, us_power_nukes_numeric),
       booktabs = T,
       caption.above = T,
       caption = "OLS Results: Sense of American Power Moderates Liberals in Iran Scenario",
       custom.coef.names = c("(Intercept)",
                             "Sense of American Power",
                             "Liberal",
                             "Male",
                             "Age",
                             "White",
                             "Republican",
                             "National Attachment",
                             "Sense of Power x Liberal"),
       stars = c(.001, .01, .05, .10),
       symbol = "\\wedge",
       label = "table:ols_iran_scenario",
       custom.model.names = c("Risky Preventive Strike",
                              "Nuclear Weapons Use"))

# --- China moderation analysis: power corrupts dovish moral foundations (MI)
china_survey$power_binary <- ifelse(china_survey$sense_power < median(china_survey$sense_power), "low", "high")
china_survey$power_binary <- factor(china_survey$power_binary, levels = c("low", "high"))

china_survey$harm_care_binary <- ifelse(china_survey$harm_care < median(china_survey$harm_care), "low", "high")
china_survey$harm_care_binary <- factor(china_survey$harm_care_binary, levels = c("low", "high"))

china_survey$fairness_recip_binary <- ifelse(china_survey$fairness_recip < median(china_survey$fairness_recip), "low", "high")
china_survey$fairness_recip_binary <- factor(china_survey$fairness_recip_binary, levels = c("low", "high"))

china_survey$individualizing_foundations <- china_survey$harm_care + china_survey$fairness_recip
china_survey$individualizing_binary <- ifelse(china_survey$individualizing_foundations < median(china_survey$individualizing_foundations), "low", "high")
china_survey$individualizing_binary <- factor(china_survey$individualizing_binary, levels = c("low", "high"))

summary(mi_china_individualizers <- lm(mi_war ~ sense_power*individualizing_binary + 
                                         ingroup_loyalty + authority_respect + nat_attach + 
                                         education + gender + age, data = china_survey))

china_survey$binding_foundations <- china_survey$ingroup_loyalty + china_survey$authority_respect
china_survey$binding_foundations_binary <- ifelse(china_survey$binding_foundations <= median(china_survey$binding_foundations), "low", "high")
china_survey$binding_foundations_binary <- factor(china_survey$binding_foundations_binary, levels = c("low", "high"))

summary(china_power_intentions <- lm(aggressive_intentions ~ sense_power*binding_foundations_binary + 
                                       nat_attach +  education + gender + age +
                                       harm_care + fairness_recip, data = china_survey))

summary(china_power_threats <- lm(us_threats ~ sense_power*binding_foundations_binary + 
                                    nat_attach +  education + gender + age +
                                    harm_care + fairness_recip, data = china_survey))

summary(china_power_intentions_table <- lm(aggressive_intentions ~ sense_power*binding_foundations + 
                                             nat_attach +  education + gender + age +
                                             harm_care + fairness_recip, data = china_survey))

summary(china_power_threats_table <- lm(us_threats ~ sense_power*binding_foundations + 
                                          nat_attach +  education + gender + age +
                                          harm_care + fairness_recip, data = china_survey))
# --- Table A10
screenreg(list(china_power_intentions_table, china_power_threats_table),
       booktabs = T,
       caption.above = T,
       caption = "OLS Results: Sense of Chinese Power Moderates Dovish Foundations in South China Sea Scenario",
       custom.coef.names = c("(Intercept)",
                             "Sense of Chinese Power",
                             "Binding Foundations",
                             "National Attachment",
                             "No College Degree",
                             "Male",
                             "Age",
                             "Harm/Care",
                             "Fairness/Reciprocity",
                             "Sense of Power x Binding Foundations"),
       stars = c(.001, .01, .05, .10),
       symbol = "\\wedge",
       label = "table:ols_china_moderation",
       custom.model.names = c("Countries Harbor Aggressive Intentions",
                              "U.S. Threatens China's Security"))

china_power_intentions_predict <- ggpredict(china_power_intentions, c("sense_power", "binding_foundations_binary")) 
china_power_threats_predict <- ggpredict(china_power_threats, c("sense_power", "binding_foundations_binary")) 

individualizers_predict <- ggpredict(mi_china_individualizers, c("sense_power", "individualizing_binary")) 

plot_df_china_individualizing_mi <- 
  rbind(data.frame(mil_int = individualizers_predict$predicted,
                   iv = individualizers_predict$group,
                   power = individualizers_predict$x,
                   ci_high = individualizers_predict$conf.high,
                   ci_low = individualizers_predict$conf.low))

plot_df_china_individualizing_mi$iv <- factor(plot_df_china_individualizing_mi$iv, levels = c("low", "high")) 
levels(plot_df_china_individualizing_mi$iv ) <- c("Low", "High")
plot_df_china_individualizing_mi$country <- "Chinese Public"
plot_china_individualizing_mi <- 
  ggplot(plot_df_china_individualizing_mi, aes(x = power, y = mil_int, group = iv)) +
  geom_ribbon(aes(ymax = ci_high, ymin = ci_low), alpha = .45, fill = "gray70") +
  geom_line(aes(color = iv), size = 2) +
  labs(y = "Militant Internationalism", x = "Sense of Chinese Power", colour = "Committment to\nIndividualizing\nFoundations") +
  scale_color_manual(values = c("#99b0b0", "#54545d")) +
  theme_bw() +
  facet_wrap(~country) +
  guides( colour = guide_legend(reverse = T)) +
  theme(axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14),
        axis.title.y = element_text(size = 16, margin = margin(r=10)),
        axis.title.x = element_text(size = 16, margin = margin(t=10)),
        legend.text = element_text(size = 13),
        legend.title = element_text(size = 13),
        strip.background = element_rect(fill = "#54544c", color = "black", size = 1.5),
        panel.background = element_rect(fill = "#f6f6f5"),
        panel.spacing = unit(1.3, "lines"),
        strip.text = element_text(size = 14, color = "gray94"),
        plot.margin = margin(23,23,23,23))

# --- Figure 3:
ggpubr::ggarrange(plot_us_ideology_mi, plot_china_individualizing_mi, nrow = 2, labels="AUTO",  font.label = list(size = 20), heights = c(.95,.98))


# --- Table A9:
# The above moral foundations results are dichotomized for ease of presentation in the main text;
# these appendix models represent MF using the original scaled values
summary(mi_china_harm_simple <- lm(mi_war ~ sense_power*harm_care, data = china_survey))
summary(mi_china_harm <- lm(mi_war ~ sense_power*harm_care + nat_attach +  education + gender + age +
                              harm_care + fairness_recip + ingroup_loyalty + authority_respect, data = china_survey))

summary(mi_china_fairness_simple <- lm(mi_war ~ sense_power*fairness_recip, data = china_survey))
summary(mi_china_fairness <- lm(mi_war ~ sense_power*fairness_recip + nat_attach +  education + gender + age +
                                  harm_care + fairness_recip + ingroup_loyalty + authority_respect, data = china_survey))

screenreg(list(mi_china_harm_simple, mi_china_harm, mi_china_fairness_simple, mi_china_fairness),
       custom.coef.names = c("(Intercept)",
                             "Sense of State Power",
                             "Harm/Care",
                             "Sense of State Power $\\times$ Harm/Care",
                             "National Attachment",
                             "Less than College",
                             "Male",
                             "Age",
                             "Fairness/Reciprocity",
                             "Ingroup/Loyalty",
                             "Authority/Respect",
                             "Sense of State Power $\\times$ Fairness/Reciprocity"),
       booktabs = T,
       custom.model.names = c("(1)",
                              "(2)",
                              "(3)",
                              "(4)"),
       caption.above = T,
       caption = "OLS Results: The Sense of Power Corrupts Dovish Moral Foundations")


# --- Experimental Results: Sense of Power Activates Militant Internationalism --- #
# This experiment randomly varies relative state power (rise, decline, no power info), finding that rise activates hawkishness (MI)
# relative to decline; and, decline activates dovishness relative to the baseline, no power info condition.

prolific_survey$power_trt <- factor(prolific_survey$power_trt, levels = c("low_power", "control", "high_power"))
summary(rise_vs_decline <- lm(mi_factor ~ power_trt + age + gender + ethnicity + education + ideology + pid + nat_attach, subset(prolific_survey, power_trt %in% c("high_power", "low_power"))))

prolific_survey$power_trt <- factor(prolific_survey$power_trt, levels = c( "high_power", "control", "low_power"))
summary(decline_vs_baseline <- lm(mi_factor ~ power_trt + age + gender + ethnicity + education + ideology + pid + nat_attach, subset(prolific_survey, power_trt %in% c("control", "low_power"))))

# --- Table A15
screenreg(list(rise_vs_decline, decline_vs_baseline),
       custom.coef.names = c("(Intercept)", "Rise Condition", "Age", "Male", "White", "Bachelors Degree", "Conservative", "Republican", "National Attachment", "Decline Condition"),
       caption = "Rise/Decline Experiment: OLS Results",
       caption.above = T,
       booktabs = T)

# plot it
mod_list <- list("rise_decline" = rise_vs_decline,
                 "decline_control" = decline_vs_baseline)
results_df <- data.frame()
for(i in 1:length(mod_list)){
  mod_i <- mod_list[[i]]
  results_df <- rbind(results_df,
                      data.frame(ci_95_high = confint(mod_i, level = .95)[,2],
                                 ci_95_low = confint(mod_i, level = .95)[,1],
                                 ci_90_high = confint(mod_i, level = .90)[,2],
                                 ci_90_low = confint(mod_i, level = .90)[,1],
                                 est = coef(mod_i),
                                 variable = names(coef(mod_i)),
                                 model = names(mod_list)[i]))
  
}
results_df$variable <- revalue(results_df$variable, c("power_trthigh_power" = "Rising US Power\n(vs. Declining US Power)",
                                                      "power_trtlow_power" = " Declining US Power\n(vs. Baseline Condition)",
                                                      "age" = "Age",
                                                      "gendermale" = "Male",
                                                      "ideology" = "Conservative",
                                                      "nat_attach" = "National Attachment",
                                                      "educationbachelors" = "College Educated")) 

results_df$model <- factor(results_df$model, levels = c("rise_decline", "decline_control"))
results_df <- results_df[!results_df$variable%in%c("(Intercept)", "ethnicitywhite", "pid"),]
results_df$variable <- 
  factor(results_df$variable, rev(c("Rising US Power\n(vs. Declining US Power)", " Declining US Power\n(vs. Baseline Condition)",
                                    "Male", "Age", "Conservative", "College Educated", "National Attachment")))
plot_col <- ifelse(results_df$variable %in% c("Rising US Power\n(vs. Declining US Power)", " Declining US Power\n(vs. Baseline Condition)"), "#54545d", "#99998e")
# --- Figure 4:
ggplot(results_df, aes(x = est, y = variable, color = plot_col)) +
  geom_vline(xintercept = 0, linetype = "solid", color = "black", size = .8) +
  geom_linerange(aes(xmin = ci_95_low, xmax = ci_95_high), size = 1, color = "black") +
  geom_linerange(aes(xmin = ci_90_low, xmax = ci_90_high), size = 2.8) +
  geom_point(size = 5, color = "black") +
  geom_point(size = 4) +
  facet_wrap(model~., scales = "free_y", nrow = 2) +
  theme_bw() +
  coord_cartesian(xlim = c(-.25,.50)) +
  labs(x = "Effect on\nMilitant Internationalism", y = NULL) +
  scale_color_manual(values = plot_col) + 
  theme(axis.title.x = element_text(size = 17, margin = margin(t=15)),
        axis.text.y = element_text(size = 16),
        axis.text.x = element_text(size = 15),
        strip.text = element_blank(),
        legend.position = "none",
        panel.spacing = unit(3, "lines"),
        panel.background = element_rect(fill = "#f6f6f5"),
        strip.background = element_blank(),
        plot.margin = margin(r=15, l=15, t= 15, b = 15))

# --- Mediation analyses 
# So, there's a main effect of treatment on MI. here, we want to compare the mechanistic pathways of
# trt --> sense of power --> hawkishness, as this paper expects, versus an alternative pathway of 
# trt --> hawkishness -->  motivated desires to view state as powerful. basically, rather than sense of power
# activating hawkishness, perhaps hawks are motivated to evaluate their state as powerful

# power --> sense power --> hawkishness pathway (i.e., the paper's primary expectation)
summary(fit_med <-lm(sense_power_factor ~ power_trt + age + gender + ethnicity + education + ideology + pid + nat_attach, prolific_survey))
summary(fit_dv <- lm(mi_factor ~ power_trt + sense_power_factor + age + gender + ethnicity + education + ideology + pid + nat_attach, prolific_survey))
med_power <- mediation::mediate(fit_med, fit_dv, treat="power_trt", mediator='sense_power_factor', boot=T, control.value = "low_power", treat.value = "high_power")
summary(med_power) # 64.4% proportion mediated

# power --> hawkishness --> sense power pathway (i.e., an alternative motivated bias explanation)
summary(fit_med <-lm(mi_factor ~ power_trt + age + gender + ethnicity + education + ideology + pid + nat_attach, prolific_survey))
summary(fit_dv <- lm(sense_power_factor ~ power_trt + mi_factor + age + gender + ethnicity + education + ideology + pid + nat_attach, prolific_survey))
med_mi <- mediation::mediate(fit_med, fit_dv, treat="power_trt", mediator='mi_factor', boot=T, control.value = "low_power", treat.value = "high_power")
summary(med_mi) # versus 5.5% proportion mediated

df_med_results <- 
  rbind(
    data.frame(
      est = c(med_power$d0, med_power$z0, med_power$tau.coef, med_power$n0),
      ci_low = c(med_power$d0.ci[1], med_power$z0.ci[1], med_power$tau.ci[1], med_power$n0.ci[1]),
      ci_high = c(med_power$d0.ci[2], med_power$z0.ci[2], med_power$tau.ci[2], med_power$n0.ci[2]),
      type = c("ACME", "Direct", "Total", "Prop"),
      condition = "rise",
      pathway = "sense power --> mi",
      stringsAsFactors=FALSE),
    data.frame(
      est = c(med_mi$d0, med_mi$z0, med_mi$tau.coef, med_mi$n0),
      ci_low = c(med_mi$d0.ci[1], med_mi$z0.ci[1], med_mi$tau.ci[1], med_mi$n0.ci[1]),
      ci_high = c(med_mi$d0.ci[2], med_mi$z0.ci[2], med_mi$tau.ci[2], med_mi$n0.ci[2]),
      type = c("ACME", "Direct", "Total", "Prop"),
      condition = "rise",
      pathway = "mi --> sense power",
      stringsAsFactors=FALSE)
  )
df_med_results$pathway <- factor(df_med_results$pathway,
                                 levels = c("sense power --> mi", "mi --> sense power"),
                                 labels = c(expression("Sense of Power" %->% "MI"),
                                            expression("MI" %->% "Sense of Power")))
# --- Figure A2
ggplot(subset(df_med_results, type != "Prop"), aes(x=est, y = condition, group = type)) +
  geom_vline(xintercept = 0, linetype="solid", color = "gray10", size=.8) +
  geom_linerange(aes(xmin = ci_low, xmax = ci_high), size = 1.7,  position=position_dodge(width=-0.55), color = "gray20") +
  geom_point(aes(fill = type), pch = 21, size = 7, position = position_dodge(width = -0.55)) +
  theme_bw() +
  facet_wrap(~pathway, labeller = label_parsed) +
  labs(x = "Estimate") +
  scale_y_discrete(labels = c("Rise Condition\n(vs. Decline Condition)")) +
  scale_fill_manual(values=c("#99998e", "#d4cac9", "gray26")) +
  theme(axis.text.y=element_text(size=20, color="black"),
        axis.text.x=element_text(size=20),
        legend.title=element_blank(),
        legend.position="top",
        legend.text = element_text(size=21),
        axis.title.y = element_blank(),
        axis.title.x = element_text(size = 22, margin = margin(t = 14, r = 10, b = 0, l = 0)),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(linetype = "dashed", color = "gray10", size = .5),
        plot.title = element_text(hjust = 0.5),
        plot.margin = margin(r=10, l=10, b = 10, t = 10),
        strip.text = element_text(size = 20),
        panel.spacing = unit(1.5, "lines"),
        panel.background = element_rect(fill = "#f6f6f5"),
        strip.background = element_rect(fill = "gray90", color = "black", size = 1.5))

